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> \ Abstract 
m 

' We address the issue of mobility of localized modes in two-dimensional nonlinear Schrodinger 



lattices with saturable nonlinearity. This describes e.g. discrete spatial solitons in a tight-binding 
approximation of two-dimensional optical waveguide arrays made from photorefractive crystals. 
We discuss numerically obtained exact stationary solutions and their stability, focussing on three 



C ' different solution families with peaks at one, two, and four neighboring sites, respectively. When 

> 

varying the power, there is a repeated exchange of stability between these three solutions, with 



^ \ symmetry-broken families of connecting intermediate stationary solutions appearing at the bifur- 

cation points. When the nonlinearity parameter is not too large, we observe good mobility, and 
a well defined Peierls-Nabarro barrier measuring the minimum energy necessary for rendering a 
stable stationary solution mobile. 
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I. INTRODUCTION 



There is a large current interest in nonlinear mechanisms for storage and transport of 
localized coherent packages of energy in spatially periodic systems, which may be described 
by discrete lattice models in one, two, or three spatial dimensions. Such systems quite gener- 
ally may support intrinsic localized modes, or discrete breathers, which are exact, generally 
time-periodic, spatially localized solutions to the underlying dynamical lattice equations 
and, under certain conditions, also may show mobility. For a recent review of properties and 
applications of such intrinsic localized modes, see In particular, within nonlinear optics 
such localized modes may appear as discrete spatial solitons in periodic structures describing 



arrays of coupled waveguides. See, e.g. 



for reviews of experimental observations as 
well as theoretical modelling of discrete spatial solitons, and J,|5| and references therein for 
more recent experimental results. 

In such systems, the longitudinal variable along the waveguides plays the role of time in 
the dynamical lattice equations, while the arrays themselves may be either one-dimensional 



:iD) (e.g. 




or two-dimensional (2D) (e.g. ^) transversally. 
Traditionally j6|, |l^ , most attention has been put on waveguide arrays constructed from 
Kerr nonlinear media, in which case the appropriate lattice model derived from coupled- 
mode theory (tight-binding type approximation) is the cubic discrete nonlinear Schrodinger 



(DNLS) equation (see, e.g., |ll 
However, more recently (e.g. 



12| . In particular, optically-induced lattices with 



for a review on properties and applications of this model), 
y, y, y, Ol and references therein) much experimental effort 
has been devoted to generating discrete solitons in photorefractive media, where the nonlin- 
earity is not anymore of pure Kerr type 

focusing as well as defocusing nonlinearities were created e.g. in [7|, |8|,|9[, leading to direct 
observation of many earlier predicted phenomena such as spatial gap solitons [3] and solu- 
tions of different symmetries (odd and even) in ID arrays [n addition, the recent works 

m 

lave reported observation of spatial gap solitons and self-trapping also in ID perma- 
nent waveguide arrays with photovoltaic defocusing nonlinearity, showing strong saturation 
at higher powers. In this case, the appropriate lattice model in the tight-binding limit is a 
ID DNLS equation with saturable on-site nonlinearity (a discrete version of the Vinetskii- 

17, l3- 



15 



16 



Kukhtarev [13| equation), which was introduced in 14] and further studied in 
It is important to remark, that the discrete model used in ^| to model waveguide arrays 
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with photovoltaic defocusin g n onhnearity is mathematically equiyalent to the model for 
focusing nonlinearity of |l4j|l5|- A recent discussion and experimental demonstration of this 
equiyalence was eiven in ^1. 

One remarkable property of the ID saturable DNLS equation, discovered in [15||, is the 
boundedness, and at certain points even vanishing, of the so-called Peierls-Nabarro (PN) 
potential barrier as a function of the soliton power. The PN barrier for DNLS-like equations 
may be defined as the difference in energy (Hamiltonian) between the two fundamental 
localized modes centered at respectively in-between lattice sites, at the same power (norm) 
(the latter being a conserved quantity for the dynamics). Thus, it is expected to give a 
lower bound to the amount of additional energy necessary to render a stable stationary so- 
lution mobile. It was also numerically confirmed 17, 1^ that the mobility of high-power 
localized solutions in the saturable ID DNLS model was considerably enhanced compared 
to the cubic DNLS model, where the PN-barrier grows rapidly as the power increases js^. 
Another peculiar property of the ID saturable DNLS model is the existence of a family 
of exact analytical sech-shaped solutions, centered at arbitrary lattice positions, for some 
particular regimes of parameter values The power of these solutions, which are lin- 

early marginally stable, depends continuously on the position of their center, and it can be 
shown that this solution family exists in a neighborhood of the first zero of the PN barrier, 
bifurcating with the site-centered and bond-centered solution families exactly at the points 
where these families exchange their stability. Thus, the analytical solutions constitute a fam- 
ily of 'intermediate' symmetry-broken solutions connecting site-centered and bond-centered 
solutions, analogously to the scenario for enhanced mobility described for another type of 



extended ID DNLS equation in 



m 
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21[ , as well as for more general chains of coupled oscillators 



However, for 2D arrays, much less is known about discrete soliton mobility. Essentially, as 
was described in [2^, for cubic nonlinearity only wide discrete solitons are mobile. However, 
these are unstable to a 'quasicollapse' process, so that after moving a few lattice sites the 
broad solitons self- focus into narrow localized peaks, which get pinned by the lattice 3]. It 
is our purpose here, to show that a saturable nonlinearity may, under certain conditions, lead 
to highly mobile discrete solitons also for 2D waveguide arrays. The study of discrete soliton 
mobility has been suggested to be one of the most important issues in the implementation 
of these theoretical concepts in all-optical switching schemes for one-dimensional nonlinear 
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arrays and very recently for two-dimensional ones in a 'reduced' geometry [25 1. For 
this reason, the two-dimensionality can be viewed as a large improvement in this direction, 
because of the promising possibility of multiple-site connections and the direct integration 
with photonic crystals. 

An interesting analogy may be drawn to what is known within the field of polarons, where 
in the case of the standard semiclassical Holstein model with harmonic on-site oscillator 
potentials the stable (small) polarons are always pinned to the lattice in 2D, while if a 
realistic saturable anharmonicity is taken into account for the oscillator potentials, moving 

n 

polarons may also exist .2^ . In the former case, the static polarons are obtained as stationary 
solutions to the cubic DNLS equation, while in the latter case as solutions to the saturable 
DNLS equation (although the dynamics of the polarons is more complex as it involves 
electron as well as phonon degrees of freedom). 

The structure of this paper is as follows. In Sec. Ullwe describe the 2D saturable DNLS 
equation, its basic properties, and the region of existence for localized solutions. Next, in Sec. 
mil we analyze the stability of these solutions and we introduce the concept of Intermediate 
Solutions. In Sec. llVt we present our results for the mobility of localized solutions in 2D 
arrays. Finally, Sec. IVl concludes the paper. 

II. MODEL 



We consider the following (general) form of 
isotropic medium, analogous to the ID model in 



;he 2D saturable DNLS equation for an 



15 



3, 



'^ + ^""--" (l + |«„.„P) =°- 

where ^ is the normalized propagation distance, Un,m describes the electric field amplitude 
in the {n,m} site, A represents the 2D discrete Laplacian, Aun,m = ""n+i.m + Un-i,m + 
Un,m+i + Un,m-i- The parameter 7 is given by the ratio between the nonlinear parameter 
and the coupling constant 14 , Q] . We choose 7 > in our computations without loss 
of generality. Note that, although this implies that we are formally restricting to a focusing 
nonlinearity, our results are immediately translated to the defocusing case through the stag- 
gering transformation Un,m {—^)"^~^"'Un,m,^ We use an isotropic approximation 
which essentially considers the coupling between neighboring sites in the h and rh directions 
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as equal. In the experiment, the intrinsic anisotropy of photorefractive materials can be re- 
duced by changing the lattice orientation relative to the crystal axis j27|. The two conserved 
quantities for (^Q) are the energy (Hamiltonian) 

N,M 

H = — [(ifn+l.m + Un,m+l) U*n,Ta 



n,m=l 



■-ln{l + \Un,m\'^) + C.C. , (2) 



and the Power (norm) 



N,M 



n,m=l 

Thus, for small P, Eq. reduces to the cubic DNLS equation with focusing nonlinearity 
of strength 7, which then may be replaced by unity through rescalings, while for larger P 
saturation effects become important and 7 is left as an independent parameter. Note that 
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id, 



slightly different forms of the saturable DNLS equation were considered (in ID) in 
but these are easily shown to be equivalent to the model of [iJ, ll^, lla through simple 

n n n n 

rescalings and gauge transformations [1^. In [U, [l^ 112I| the authors study analytically 
the case of 7 ~ 9, based on experimental parameters. We decided to study, for the two- 
dimensional array, a similar case for 7 = 10. As it further will be explained, we also study 
the case of 7 = 4 where we found that a good mobility can be observed for lower level of 
power (P). 

To study the dynamics we introduce a definition of the soliton center in the direction n 
(horizontal direction) function of ^, 

N,M N,M 

< n > (0 = KAOW E l^nAOl'. (4) 

n,m=l n,m=l 

with an equivalent definition in the vertical direction m for < m > (^). We define "axial 
propagation" as the propagation in the n {k^ 0,ky = 0) or in the rh [ky 7^ 0, /Cj, = 0) 
directions, and "diagonal propagation" when we launch the solution with the same angle in 
both directions {kx = ky 0). 

Stationary solutions to ((T)) are those of the form Un,m{0 = Un,mC'^^^, where A represents 
the spatial frequency in the propagation direction. As we will mainly consider the funda- 
mental localized solutions, we will assume to be real (although stationary localized 
solutions with nontrivial complex Un,m, such as vortex solitons, generally also exist in 2D 




FIG. 1: (Color online), (a) 00 profile, (b) OE profile, (c) EE profile, (d) and (e) P versus A 
for 7 = 4 and for 7 = 10, respectively. Diamonds, stars, and squares represent 00, OE, and EE 
solutions, respectively. System size = M = 15 ((d) and (e) insets: N = M = 31). 

photorefractive optical lattices |2S|). Using a Newton- Raplison method, we have found the 
same types of stationary solutions as_described, e.g., in Ref. j^l for the 2D cubic model. 



Using the usual optical terminology ^, we rename these solutions: one-peak solutions, as 
Odd-Odd (00) (Fig. Ufa)); two-peaks solutions, as Odd-Even (OE) (Fig. [Tfb)); four-peaks 
solutions, as Even- Even (EE) (Fig. ^c)). The region of existence (see Fig. ^) of localized 
solutions in this system has a very different structure compared to the cubic DNLS case 
This can be understood by considering the properties of plane wave-solutions to 
Eq. (PJ) in the limits of small and large amplitudes, respectively. For small-amplitude plane 
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waves, it is easy to show, by neglecting the term [Mn^mP in the denominator of the last 
term in Eq. (P), that the linear band corresponds to {—7 — 4, —7 + 4}; the superior limit 
(A = for 7 = 4 and A = —6 for 7 = 10, see Fig. corresponding to constant-amplitude 
solutions (zero wave vector) and being the border where small-amplitude delocalized and 
localized solutions are connected. On the other hand, in the high-amplitude limit, we may 
completely neglect the last term in Eq. due to the saturable nature of the nonlinearity, 
and thus in this limit plane waves constitute the band —4 < A < 4 (independent of 7); again 
with the superior limit A = 4 corresponding to the constant-amplitude (zero wave vector) 
solution. The superior limit for localized solutions is observed to be A = 4 by increasing the 
frequency (see Fig. Q), and thus the region of existence for localized modes is between the 
low-amplitude and high-amplitude limits for the upper band edge (zero wave vector) plane 
wave. 

Note that, similarly to the cubic 2D DNLS equation [23, Q, Q Q, the power for 
the localized solutions is non-zero in the small-amplitude limit, and, when the frequency 
increases, the power first decreases until it reaches a minimum. When the frequency of 
the localized solutions increases further, the power increases indefinitely. However, the 
nonlinearity is saturable, which means that there is some threshold value for the amplitudes 
beyond which localization effects should diminish. For instance, if we take an 00 solution 
(Fig. ^a)), we can see that when the central site achieves such an amplitude threshold 
value, the surrounding sites begin to increase their amplitudes, delocalizing the profile and 
evolving to a high-amplitude plane wave. An analogous scenario was discussed, in a polaron 
context, in Ref. 12611 . and similar properties can also be described for the ID saturable 
problem jl^ . 




III. STABILITY ANALYSIS OF TWO-DIMENSIONAL LOCALIZED SOLU- 
TIONS 

Linear stability of stationary solutions may be investigated in a standard way (see, e.g., 
Ref. 0]) by writing M„,m(0 = (""n,™ + 4>n,m{0)^^^^y leading to the hnearized equation for 
(f)n,m, the perturbation function (PF). To solve this problem numerically, we use the technique 
outlined in Ref. for the stability analysis of localized solutions in a similar nonlinear 
discrete medium. We split the PF into real and imaginary parts, 0„_m = Xn,m + iyn,m 
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{x,y G 3?), and insert them into the hnearized equation. Taking a two-dimensional square 
discrete array (A^ = M) to compute stationary solutions and to perform their stability 
analysis, we proceed to map the 2D problem to a ID representation. Defining Z as the 
vector Zjq2 [Z = X + iY), we can write the equations for the perturbation functions 

in a one- dimensional representation as 

X + AY = and F - BX = ^ 

X + ABX = and F + BAF = 0, (5) 

where A and B are A^^ x A^^ matrices, depending of the Un,m profiles. Now, the stability 
analysis consists in finding the eigenvalues of the matrix AB (the matrix BA has the same 
eigenvalues If all the eigenvalues are positive the solutions of the problem (PF's) are 

oscillatory functions, which implies stability. On the other hand, a negative or complex 
eigenvalue means the existence of exponentially increasing PF's, implying instability. 

For a given frequency A we compute the Power P, the Hamiltonian H, and the AB 
eigenvalues for the 00, OE, and EE solutions. 

We take the most negative eigenvalue of each solution, calling it "G", and we plot it 
for different powers in Fig. Efa) and (b) for 7 = 4, and in Fig. ISfa) for 7 = 10. G = 
or G < implies stability or instability, respectively (there are no complex eigenvalues for 
these solutions). An oscillatory behavior in the stability of solutions can clearly be viewed 
in these figures. Depending on the level of Power, a change in stability for different solutions 
may 
Ref. 

in properties since it has three main localized solutions, and because, as we will show below, 
the exchange of stability produces good mobility of very localized solutions for 2D arrays. 
One way to understand the concept of exchange of stability between two different so- 



De observed. This result is in direct concordance with the previous result showed in 
15| for the one-dimensional problem. However, the two-dimensional problem is richer 
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lutions is to look for Intermediate Solutions (IS), as was done for ID lattices in j2l| 
Such solutions, having symmetry-broken profiles interpolating between the two solutions 
which exchange their stability, typically only exist in a limited parameter regime, connect- 
ing through bifurcations with the two solutions of higher symmetry. Generally, as discussed 
in [2^, the IS may be either linearly stable, in which case the two symmetric solutions are 
simultaneously unstable, or the IS may be unstable, when the two symmetric solutions are 
simultaneously stable. Here, we find IS connecting two solutions that share stability, and 
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c) 



I 3 3 7 8 9 II 13 15 




FIG. 2: (Color online). Results for 7 = 4. (a) Smallest eigenvalue G versus Power, (b) Zoom of 
(a) in the region P ~ 9.4. (c) Profile of intermediate solution (IS). Diamonds, stars, squares, and 
triangles represent 00, OE, EE, and IS solutions, respectively. 

thus these IS are unstable. 

For 7 = 4, we found that the first "exchange region" (defined as the region where two 
solutions are stable simultaneously) is between the 00 and the OE solutions. This region 
can be observed in Fig.|21^b), where we also present the profile for the IS (Fig.|2fc)). For this 
value of 7, we can see in Fig. |21^a) that the exchange regions when increasing the power are, 
consecutively: 00-OE, OE-EE, EE-OE, OE-00. Multiple crossing points between the three 
different solutions can be observed by plotting the Hamiltonian versus Power, coinciding 
with the exchange regions defined before. These energy crossings imply that the new stable 
solution has lower energy compared to the others. This can be interpreted as an oscillation in 
the PN barrier However, there is no formal definition of this barrier for 2D arrays, and 
it generally must depend on the direction of motion. We may define it loosely as the largest 
difference in energy, at constant norm, between two stationary solutions of the system, close 
to which a localized mode must pass when moving adiabatically through the lattice in a 
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certain direction. (Note that, even in ID, the simple definition of the PN barrier used e.ff. 

n 

in Ref. [15|] as the difference between site-centered and bond-centered solutions is generally 
not appropriate when these are simultaneously stable, since the energy of the IS may be 
considerably larger). The behavior of this barrier in arrays with saturable nonlinearity is 
remarkable. As seen from Fig. |21^a), for each of the three fundamental localized solutions, 
we can pass from a region where the solution is unstable to one where it is stable, by tuning 
the corresponding solution power (analogously to the ID scenario shown in Ref. 15[)- As 
we show explicitly in Sec. IIVI this exchange of stability properties for different solutions 
implies an improvement in the mobility for discrete arrays, which is an important aspect 
for the implementation of this concept in all-optical switching systems. (We do not here 
explicitly show the results for H versus P for 7 = 4, since the relative differences between 
the respective energies are too small to be clearly represented graphically. See the discussion 
in Sec. HVlfor some typical values.) 

Compared to the previous case, for 7 = 10 (implying a higher nonlinearity or a lower 
coupling between sites) solutions are more localized as expected. The stability diagram is 
presented in Fig. IHfa). As in the former case, multiple exchange regions where two solutions 
are stable at the same level of power can be observed: 00-OE, OE-EE, EE-00, EE-OE. It 
is very clear from Fig. Efa) that the IS (Fig. El^b)) exists in the exchange region where the 
00 and OE solutions are stable (region enclosed by a circle in Fig. Efa)). In this case the 
differences in energies are bigger, and we show in Fig. IHfc) a plot for the Hamiltonian versus 
Power. In this figure the way in which the IS connect the OE and the 00 solutions can be 
very well observed. Before the crossing point, the energy difference A/J = Hqq — Hqe is 
lower than zero. After the crossing point this difference has changed its sign. This is a clear 
confirmation that the oscillation in the PN potential, if defined in the "naive" way as AH, is 
due to the exchange in the stability of different solutions. The OE solution becomes stable 
(see Fig. Efa)), then it shares stability with the 00 solution, and finally it continues to be 
stable while the 00 solution becomes unstable. The same behavior is described by energy 
considerations from Fig. Efc). Note however, that the actual energy barrier to overcome 
for a localized solution moving adiabatically in an axial direction is larger than AH in the 
regime of simultaneous 00 and OE stability, if the path goes via the IS whose energy is 
larger. This will be illustrated in Sec. IIVI 
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FIG. 3: (Color online). Results for 7 = 10. (a) Smallest eigenvalue G versus Power, (b) Profile of 
intermediate solution (IS), (c) Energy H versus Power for the region enclosed with a circle in (a). 
Diamonds, stars, squares, and triangles represent 00, OE, EE, and IS solutions, respectively. 

IV. MOBILITY OF TWO-DIMENSIONAL LOCALIZED SOLUTIONS 

We study the mobility of 2D discrete solitons by solving numerically Eq. for initial 
conditions obtained by slightly perturbing the stationary solutions described above. To be 
specific, we here choose the 00 solution and study mobility in regions where this profile is 
always stable. In our simulations, we always consider perturbations obtained by "kicking" 
the initial 00 solutions using: Un,m(0) = Un,mGxp[i {kxii + kym)]. Note that such perturba- 
tions do not change the power, but generally increase the energy compared to the stationary 
solutions. 

In Fig. m we show dynamical results for 7 = 4 and 7 = 10. First, for 7 = 4, we 
study the dynamics for small-power solutions, P A. Fig. Efa) shows the average axial 
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position (jU for an 00 profile (belonging to the branch with dP/dX > in Fig. dd)) 
kicked in the axial direction. In this region of power the 00 mode is stable, while the OE 
and EE modes are both unstable (see Fig. I21^a)). The additional energy due to the kick 
corresponds to AH = 0.0038, which is larger than the energy difference between the 00 
and OE configurations, AH = 0.0032. For this reason, the solution first begins to move, 
but is then trapped 4 sites away from the input center in a (stable) 00 configuration. Due 
to radiation losses, it has not sufficient energy to overcome the next barrier, but instead 
begins to oscillate around its new center due to the excitation of its stable internal mode. 
An interesting detail in this result is the structure of the curve. When the 00 mode 
propagates, it changes its form between the 00 and OE profiles. It is clear from the figure 
that the maximum propagation velocity corresponds to the 00 configuration (bigger slope, 
potential well), and the minimum velocity (slope ~ 0, saddle point) corresponds to the OE 
configuration (centered between two sites in the propagation direction). Similar results for 
ID systems have been described, e.g., in Refs. |22L l33l|. 

For the same parameter values, we illustrate in Fig.|31^c) the propagation of an 00 profile 
kicked in the diagonal direction. This plot describes the identical evolution in both axes. 
It can be observed that the 00 profile propagates in the diagonal direction, changing its 
form from an 00 configuration to an EE one (see the inset). The final switching here was 
two sites in both directions, from the input position {8, 8} to {10, 10}. Here, the maximum 
propagation velocities correspond to the 00 configurations, and the almost zero velocities 
correspond to the EE configurations. As before, the soliton will not be able to continue its 
propagation because of the radiation losses. In this case, the diagonal kick corresponds to 
AH = 0.007, while the energy difference between 00 and EE configurations at this power 
is AH = 0.006. 

The dynamics shown in Figs. |^ (a) and (c) corresponds essentially to that of a typical 
ID cubic DNLS case in the regime of low power, where the PN barrier is relatively small, 
and the site-centered solution possesses a symmetry-breaking internal translational mode 
(see, e.g., js^), which may be excited by kicking the solution. Now, if we compare with the 
2D cubic DNLS case, the dynamics and the regions of existence of the stationary localized 
solutions are very different. First, in a P versus A diagram for the 2D cubic DNLS case [2^, 
there are different power thresholds for 00, OE, and EE solutions, far away separated in 
power. This implies that below a certain level of the power, it is not possible to move the 
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00 profile since no OE or EE solution exists at the same power. If we consider a higher level 
of the power (on the high-frequency branch with a stable 00 solution), where two or three 
solutions exist, the profiles are too localized and do not possess any localized symmetry- 
breaking internal modes. This implies also that the PN potential is large, and mobility is 
not possible by kicking the solutions since there are no internal translational 'depinning' 
modes to excite. Thus, for cubic 2D DNLS, mobility is possible only on the 'continuumlike' 
low-frequency branch with dP/dX < 0, where 00, OE, and EE solutions are all unstable, 
but such moving solutions will typically 'quasicollapse' into localized pinned solutions as 
described in |23]. 

In the 2D saturable case, these thresholds still exist but are much closer in power. For 
7 = 4, the power thresholds for different solutions are: Pqo ~ Pqe ~ -Pee ~ 3.07 (cf. 
Fig. ^d)). This means that all solutions have essentially the same power threshold value, 
and that we are in principle able to observe good mobility in the low-power regime of the 
branch with dP/dX > for P > 3.08. For 7 = 10, the differences in power are higher than 
in the previous case, but they are still rather small. The corresponding threshold values 
are (see Fig. life)): Pqo = 0.90, Pqe = 1-13, and Pee = 1-17. This behavior is another 
remarkable property of the saturable nonlinearity. As is wellknown (see, e.g., js^), in two- 
dimensional discrete nonlinear systems with an effectively cubic (or stronger) nonlinearity, 
a power threshold value always exists for localized solutions. Our results suggest that, 
increasing the value of 7, it is possible to decrease this power threshold value towards zero 
(for 7 = 100 this value is Pqo = 0.059). In the small-amplitude regime, the dynamics is 
essentially governed by the cubic term in the Taylor expansion of the saturable nonlinearity, 
and thus Eq. (P) is well approximated by the cubic DNLS model with effective nonlinearity 
parameter 7P. Thus, the threshold power should scale as 7"^ for large 7. A similar result 
was also discussed in Ref. j3| (see, e.g.. Fig. 6 in this paper). 

Now, we go further to study the exchange regions where we observe "multistability" of 
solutions. In the region of power shown in Fig. Efb), we study the first exchange region 
between the 00 and the OE solutions. In Fig. IHd) we show the effect of kicking the 
00 solution in the axial direction. The energy added to the profile due to the kick is 
AH = 0.0003, and the energy difference between the 00 and the OE configurations is 
AH = —0.0001, i.e., the OE solution has the lowest energy. However, in this case there 
exists also an intermediate solution (Fig. Efc)), which is important to explain the mobility 
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we observe in this dynamics. The energy difference between the 00 and the IS solutions is 
/\H = 0.0002. In fact, the solution can move across the array as far as the radiation losses 
permits it. The regions in which the profile changes its velocity can clearly be observed 
in this figure. First, for the 00 and the OE solutions the velocity has maxima, which is 
corresponding with the stability analysis where both solutions are stable. Both solutions 
correspond to a potential well in a dynamical representation. We can also note, that the 
velocity is larger at half-integer values than at integer values, consistent with the fact that 
the OE solution has lower energy than 00. The minimum velocities (slope ~ 0) clearly 
correspond to the IS between the 00 and the OE profiles j^^. For this case IS have shown 
to increase mobility, essentially because all solutions are very close in energy and just a 
little kick, corresponding to approximately the energy difference between the 00 and IS, is 
required to move them. An analogous scenario would be seen by using the kicked stable OE 
solution as initial condition. 

Finally, we study the case for 7 = 10. We first look for mobility in the same region of 
power as in the previous case (P ^ 4 and P ^ 9). For these powers we were not able to 
find mobility of localized solutions, essentially because the stationary solutions are distant 
in energy values. 

Then, we study the first crossing point between the energy of the 00 and the OE solutions 
(Fig. Etc)). From the previous discussion, it is clear that in the exchange regions, the 
mobility of solutions depends on the IS. So, if we want to switch the profile in the lattice, 
we first have to overcome the energy barrier between the 00 and the IS. We study two 
cases where we found some mobility. First, we take an 00 profile for a power lower than 
the crossing point power (P ^ 21.5), i.e. Pqo = 20.43. For this power both solutions are 
stable, but the 00 solution has lower energy than the OE solution. Therefore, we expect 
it to be most probable to have an 00 stable configuration at the end of the process by 
switching the profile. However, if radiation losses are large, the solution might in principle 
also be trapped in a metastable OE configuration, since it might not have enough energy to 
overcome the barrier between OE and 00 created by IS. In Fig. IHf), we show an example 
of switching of an 00 mode with one site along the axial direction. The energy due to the 
kick was AH = 2.147, the energy difference between the 00 and the OE configurations was 
AH = 0.370, and the energy difference between the 00 and the IS was AH = 0.672. The 
figure shows how the 00 solution starts to move very fast (large slope), but it is not able 
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to continue its movement in the array because of the high radiation losses of the switching 
process. This is a remarkable example of discrete switching in nonlinear lattices. 

Now, we consider the situation where the 00 power (Pqo = 24.76) is higher than the 
crossing point power. Here, both 00 and OE solutions are stable, but the OE solution has 
lower energy, so from energetic arguments we would expect the most likely final profile to be 
an OE mode (however, depending on the amount of radiation losses, either an OE or an 00 
profile may be observed as final state as discussed above). The dynamics resulting from a 
kick of the 00 mode producing a AiT = 2.045 is illustrated in Fig.|3fg). The corresponding 
energy difference between the 00 and the OE modes is AiJ = —1.120, and between the 00 
mode and the IS, AiJ = 0.118. It is clear from these numbers that we need to supply the 
00 mode with some extra energy to move it across the array even though the energy for the 
OE solution is lower. This property is due to the oscillating behavior of the Hamiltonian. 
Without such crossing points, 2D mobility for high power solutions can generally not be 
expected because of the large PN barrier. Fig. E^g) shows the switching of an 00 solution 
with two sites in the axial direction from the input position. First, the 00 solution has 
sufficient energy (big slope) to overcome the IS energy barrier. Then, its velocity increases 
as it passes the OE position, and then again it slows down when approaching the next IS 
(with central position close to site 9). Now, it does not have enough energy to overcome this 
barrier, and instead it makes one full oscillation around the OE position, limited by the two 
unstable IS. However, during this oscillation it has recovered some of its energy, and can 
now pass the IS barrier to the position of the next site, where it gets temporarily trapped 
into an 00 configuration. The radiation losses are still not too high, and the 00 mode 
being energetically unstable decays after a short distance ^ into a state of large-amplitude 
oscillations around the energetically stable configuration, the OE mode centered in {9.5, 8} 
(see inset in Fig. |H^g)). Now the solution has enough energy to oscillate between two sites (9 
and 10). This oscillation produces more radiation losses, and the final profile is an 00 mode 
centered in the site number 10, trapped by the barrier created by its nearby IS. Figs. Iljb), 
(e), and (h) show different profiles in the different regions of mobility. It is clear from these 
figures that the mobility is enhanced for less localized solutions as expected from the PN 
potential concept, but it is not forbidden for highly localized solutions as we have shown. 

We have checked these dynamical cases also for a bigger array, where = M = 21. The 
quantitative results shown in Fig.|3]change somewhat due to the extension of the system, but 
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the qualitative picture for the sohton's mobihty is preserved. We set in our computations 
periodic boundary conditions, which imphes that some radiation comes back perturbing the 
final state of the profile. Therefore, the results are dependent of the boundaries. However, 
the important issue we want to emphasize here is the mobility of highly 2D localized states, 
which is independent of the array size for saturable nonlinear media. 

V. CONCLUSIONS 

In conclusion, we have analyzed in some detail the mechanisms leading to mobility of lo- 
calized modes in a two-dimensional DNLS-type lattice with saturable nonlinearity. From a 
practical point of view, the most important result from our study is the drastic enhancement 
of the mobility resulting from the saturable nature of the nonlinearity, both for low-power 
and high-power excitations. This effect, which should be observable for two-dimensional 
waveguide arrays, is in our opinion more remarkable than the similar effect previously re- 
ported in ID [l^, since stable localized excitations for pure Kerr nonlinear media are known 
to be essentially immobile in 2D. Thus, the saturability of the nonlinearity introduces new 
possibilities for power controlled steering and switching also for 2D arrays. As we have 
shown, mobility is not restricted to axial directions, but also steering in diagonal directions 
is possible. 

In principle, the choice of a smaller value of 7, for instance 7 = 8, could also be good 
for 2D mobility, because it is expected that the first stability exchange region occurs for 
lower level of power. However, it is important to notice that a decrease in the value of 7 
(see Fig. also implies an increasing value for the power threshold of the solutions, due to 
a balance between the nonlinearity and the coupling terms. We confirmed this behaviour 
for 7 = 3, where the power threshold was found to be P ^ 4.84, while the first stability 
exchange region was found to occur at P ~ 6.98. On the contrary, if we increase further the 
value of 7, for instance 7 = 20, we expect that the powers for the first exchange region will 
be much higher and, in this sense, it is not really interesting from the optical application 
point of view, which always requires low level of power. 

From a more fundamental point of view, our results also give a deeper understanding 
for the mechanisms for mobility of localized modes. In particular, this concerns the relation 
between regimes of exchange of stability between site-centered and bond-centered stationary 
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solutions, and points of vanishing of a so-called Peierls-Nabarro (PN) barrier defined naively 
as the difference between such solutions at constant norm. As we have shown, this definition 
is not appropriate in regimes where both these solutions are simultaneously stable, due to 
the existence of unstable intermediate solutions (IS) of higher energy. However, redefining 
the PN barrier as the energy difference to the relevant IS gave good agreement with the 
numerically observed additional energy necessary for making a stationary solution mobile. 
An analogous scenario should exist in ID, and gives an intuitive ex pla nation to why, in 
spite of the repeated vanishing of the "PN potential" reported in Ref. jl5| at several critical 
powers, the mobility is good only around the first of these [2^. At the other (larger) critical 
powers, the barriers created by the IS are expected to be very large, and no smooth path in 
phase space passing simultaneously close to all three stationary solutions is likely to exist. 

It is also interesting to mention, that although, in ID, both the existence of a particular 
class of IS (which were even obtained analytically) and the "PN barrier vanishing" 
were previously known for the saturable potentials, the connection between these results, 
and their relation to exchange of stability, seems so far to have gone unnoticed. 
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FIG. 4: (Color online). Mobility results. 7 = 4, P = 3.99: (a) axial propagation for kx or 
ky = 0.0323; (b) 00 profile; (c) diagonal propagation for k^ = ky = 0.031. 7 = 4, P = 9.42: (d) 
axial propagation for kx or ky = 0.00624; (e) 00 profile. 7 = 10: (f) axial propagation for kx or 
ky = 0.6 (P = 20.43); (g) axial propagation for kx or ky = 0.5 (P = 24.76); (h) 00 profile. System 
size = M = 15, periodic boundary conditions. 
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